// SparseArray.h: interface for the SparseArray class.
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_SPARSEARRAY_H__5636894A_624C_11D1_93CA_006097737B41__INCLUDED_)
#define AFX_SPARSEARRAY_H__5636894A_624C_11D1_93CA_006097737B41__INCLUDED_

#if _MSC_VER >= 1000
#pragma once
#endif // _MSC_VER >= 1000

#include <iostream>
using namespace std;
#include <cstdlib>
#include <cmath>
#include "nr_c/nr.h"

// Class for a sparse array NOT matrix. You can form a matrix out of
// an array of SparseArray. 
class SparseArray
{
  // take the dot product of a SparseArray with a another array
  friend double operator*(const SparseArray &X, const double *y);
  friend ostream& operator<<(ostream&, const SparseArray&);
public:
	SparseArray(double i_Limit = 1e-8)	{ N=0; Limit = i_Limit; }
	virtual ~SparseArray()	{ }

	void    operator= (SparseArray &X);
	double &operator[](const int   &m);

	void    Clean();

  void multiplyElementByElement(const double *y);
  void multiplyElementByElement(const double y);

private:
	int N; // number of elements in sparse array
  int NonZero[32];  // sparse array index numbers
	double Value[32];  // sparse array values where Value[1] is the value of element NonZero[1]
	double Limit;
};


//////////////////////////////////////////////////////////////////////
// overloaded operators for SparseArray class
//////////////////////////////////////////////////////////////////////

inline void SparseArray::operator=(SparseArray &X)
{
N = X.N;

CopyA(Value,   X.Value,   N);
CopyA(NonZero, X.NonZero, N);
}


inline double operator*(const SparseArray &X, const double *y)
{
double Answer = 0.0;
 
int p;

for(int n=1;n<=X.N;n++)
	{
	p       = X.NonZero[n];
	Answer += X.Value[n]*y[p];
	}

return (Answer);
}

inline void SparseArray::multiplyElementByElement(const double *y)
{
  for(int n=1; n<=N; n++)
	{
	  int p = NonZero[n];
	  Value[n] *= y[p];
	}  
}

inline void SparseArray::multiplyElementByElement(const double y)
{
  for(int n=1; n<=N; n++)
	{
	  Value[n] *= y;
	}  
}

inline void SparseArray::Clean()
{	if(N)	if(fabs(Value[N])<=Limit)	N--;	}


#endif // !defined(AFX_SPARSEARRAY_H__5636894A_624C_11D1_93CA_006097737B41__INCLUDED_)
